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^ ■ Abstract 



We study the equilibrium properties of a flexible homopolymer where consecutive monomers 
are represented by impenetrable hard spheres that are tangent to each other, and non-consecutive 
monomers interact via a square-well potential. To this aim, we use both replica exchange canonical 
simulations and micro-canonical Wang-Landau techniques for relatively short chains, and perform 
a close comparative analysis of the corresponding results. These investigations are then further 
exploited to reproduce, at a much shorter scale and, hence, computational effort, the phase diagram 
previously studied with much longer chains. This opens up the possibility of improving the model 
and introduce specificities typical, among other examples, of protein folding. 
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I. INTRODUCTION 

Square-well (SW) potential has a long and venerable tradition in simple liquids and 
has become a paradigmatic test-bench for more sophisticate new approaches. Initially it 
was used as a minimal model, alternative to Lennard- Jones potential, in early attempts of 
molecular dynamics simulations of liquids 3], because it could be more easily implemented 
in a simulation code, and yet contained the most important features of pair potential for a 
liquid. Indeed, it displays both a gas-liquid and liquid-solid transitions in the phase diagram, 
with results often quantitatively in agreement with real atomistic fluids 3)]. For sufficiently 
short-range attraction, the gas-liquid transition becomes metastable and gets pre-empted by 
a direct gas-solid transition 4|, |5| • Several variants of the SW model have also beenproposed 
over the years in the framework of molecular fluids 6] and colloidal suspensions {7]. 

In the framework of polymer theory, the model is relatively less known, but it has experi- 
enced a re-surge of interest in last two decades as a reasonable compromise between realism 



and simplicity [8Hl2j. In this model, the polymer is formed by a sequence of consecutive 
monomers, represented by impenetrable hard-spheres, so that consecutive monomers are tan- 
gent to one another, and non-consecutive monomers additionally interact via a square-well 
interaction. The model can then be reckoned as a variation of the usual freely-jointed-chain 



131 ]. with the additional inclusion of a short-range attraction between different parts of the 
chain and excluded volume interactions. 

In spite of its simplicity, this model displays a surprisingly rich phase behavior, including 
a coil - globule and a globule - crystal transitions, that are the strict analog of the gas- 
liquid and liquid-solid transitions, respectively. Interestingly, even in this case a direct coil - 
crystal transition is found for sufficiently short range attraction, pushing this analogy with 
the above direct gas-solid freezing transition even more an 

The SW polymer model can also be easily adapted to mimic the folding of a protein, 
rather than a polymer. The crucial difference between polymers and proteins stems from 
the specificities of each amino acids forming the polypeptide chain that, along with the 
steric hindrance provided by the side chains, drastically reduces the number of possible 
configurations of the folded state [ijj]. As a result, one obtains a unique native state, rather 
than a multitude of local minima having comparable energies. The simplest way to introduce 
the selectivity defined by different amino acids is given by partitioning the monomers in two 



classes, having hydrophobic (H) and polar (P) characters. Under the action of a bad solvent 
and/or for low temperatures, the H monomers will tend to bury themselves inside the core 
of the globule, in order to prevent contact with the solvent (typically water). The HP model 



has been shown very effective in on-lattice studies |15Nl7l|. to describe the folding process, 
at least at qualitative level. The SW model is also reminiscent of Go-like models 18M21| 
routinely adopted in protein folding studies, where the amino acid specificities are enforced 
by including the native contact list into the simulation scheme. 

One of the main difficulties involved in numerical simulations of polymer chains, stems 
from the very large computational effort necessary to investigate its equilibrium proper- 
ties for sufficiently long polymers. This is true both using conventional canonical techniques 
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23J, and more recently developed micro-canonical approaches, such as the Wang-Landau 



method |24j. Even in the simple SW polymer, while high temperature behavior poses little 
difficulties, low-temperature/low-energy regions are much more problematic, and yet most 
interesting. With canonical ensemble simulations the system frequently gets trapped into 
metastable states at low temperatures, and with the Wang-Landau method the low tempera- 
ture results strongly depend on the lowest (ground) state energy definition, and an extension 
of the value of ground energy state to lower values requires increasingly large computational 
effort. 

It is then of paramount importance to investigate the possibility of using such models 
for shorter chains and to make a critical assessment on the reliability of the corresponding 
results. The present paper presents a first step in this direction. More specifically, our 
aims are two-fold. First, we will perform a parallel investigation of the SW model for 
relatively short chains (up to 32 monomers) using both replica exch ang e canonical Monte 



Carlo simulations 23] and Wang-Landau micro-canonical technique 24j. A second goal of 
the paper is to critically assess the possibility of inferring the full phase diagram, even for 
these relatively short chains. 

The remaining of the paper is organized as follows. In Section [Til we will introduce the 
model and the relevant thermodynamical quantities. Section III II will be devoted to a brief 
recall of the Monte Carlo simulation techniques, and Section ITVl to the obtained results. The 
paper will end with some conclusions and perspectives in Section [V] 



II. POLYMER MODEL AND THERMODYNAMICS 

Following the standard approach fill . I12I , we model the system as a flexible homopolymer 
chain formed by a sequence of iV monomers, located at positions {r 1; . . . , r^r}, each having 
diameter a (see Fig. [T]). Consecutive monomers are connected by a tethering potential 
keeping the iV — 1 consecutive monomers at fixed bond length /. Non-consecutive monomers 
are subject to the action of a square-well (SW) potential 



+00 , r < a 
— e, a < r < \a 
0, r > Xa 



(2.1] 



where r,- 



rA, and A — 1 is the well width in units of a. Here e defines the 



well depth and thus sets the energy scale. The model has a discrete spectrum given by 
E n = —en, where n is the number of SW overlaps that, in turn, depends upon A. In the 
present paper we will use l/a = 1 and values of A in the range [1.03, 1.6], although the case 
Ijo > 1 has proven to be interesting 26j too. In micro-canonical approach, a central role is 
played by the density of state (DOS) g{E) that is related to the micro-canonical entropy 



S (E) = k B lng (E) 



(2.2) 



(Izb is the Boltzmann constant) and hence to the whole thermodynamics. Here additional 
interesting quantities to infer the character of the transition are the inverse micro-canonical 



temperature 



23] 



HE) = (k B T(E)r = ^± 



and its derivate 



l{E) 



dp (E) d 2 S (E) 



dE dE 2 

Canonical averages can also be computed using the partition function 



Z(T) = J2g(E)e- E ^ T \ 



(2.3) 



(2.4) 



(2.5) 



and the probability function 



P(E,T) 



Z{T) 



9{E)e 



-E/(k B T) 



(2.6) 




FIG. 1. A homopolymer with square-well potential. The chain can be represented as a set of 
identical connected hard-spheres (diameter a) with bond lengths I = a. Non-consecutive spheres 
interact via a square- well potential of range Act. Consecutive spheres are tangent one another but 
cannot inter-penetrate (when l/a = 1). 

to find system in a conformation with the energy E and temperature T. Helmholtz and 
internal energies can then be obtained as 



F (T) = —ksT In Z (T) 

u(t) = (e) = J2ep(e,t) 



(2.7) 
(2.8) 



where the average (...) is over the probability (12.61) . 

As we shall see, two important probes of the properties of polymers are given by the heat 
capacity 

dU(T) (E 2 ) — (E) 2 



C(T) 



dT 



knT 2 



(2.9) 



and by the the mean squared radius of gyration 

1 N 

R l = JjJ2^~ t ^) 2 , (2-10) 

i=i 

where r cm = J2f=i r « * s the center of mass position. In practice, one constructs the 
P(Rg, E) distribution of R g at a given E to get the micro-canonical average (■ ■ -)e over this 
distribution 

and then the canonical average 

( r K t )) = Y,( r I)e p ( e ^- ( 2 - 12 ) 

E 

It is important to remark that derivation of thermodynamics from the g(E) is quite general, 
and does not depend on the specific method of simulation. This is therefore the optimal 
tool to compare different methods and assess the pros and cons of each of them, that is one 
of the aims of the present work. 



III. MONTE CARLO SIMULATIONS 



In this section, the equilibrium properties of the above model for homopolymer will be 
computed using micro-canonical and canonical Monte Carlo simulations, for chain lengths 
up to N = 32. This will allow us to ascertain some specificities of each of them, and to 
recover some known results at a much lower computational effort. 



A. Micro-canonical approach: Wang-Landau method 



Following general established computational protocols 22J, |23j, and refs. Ill, [12J for the 
specificities related to the polymers, we use Wang-Landau (WL) method (24J to sample 
polymer conformations according to micro-canonical distribution, by generating a sequence 
of chain conformations a — > b, and accepting new configuration b with the micro-canonical 
acceptance probability 



(a->6) 



mm 



w b g(E a ) 
w a g{E b ) 



(3.1) 



where w a and Wb are weight factors ensuring the microscopic reversibility of the moves. 

A sequence of chain conformations is generated using a set of Monte Carlo moves, which 
are accepted or rejected according to Eq. 13.11 At a randomly chosen site(s) we apply with 
equal probability pivot (acting on valence and torsional angles, chosen randomly), repta- 
tion, crankshaft or backbite (sometimes also referred to as end-bridging) moves. For the 
crankshaft Wb/w a = 1 always. Those moves generating conformations through valence angle 
random sampling must include the Wb/w a = sin# fe / sin# a ratio in the acceptance probability. 
This is the case, for instance, of the pivot move acting on valence angle, as well as of the 
reptation move, that includes random generation of valence angles at the ends. The backbite 



move 
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12] is a non-local move where one of the chain ends is chosen randomly and all 



the \x monomers, lying within a 2cr range, are counted. If \i = 0, the move fails, otherwise 
one of these neighbors is randomly selected, the bond preceding is broken, and a new bond 
is formed between the chosen neighbor site and end. As this bond length is usually larger 
than a, an additional shift is often required to satisfy the fixed bond constraint, and the 
re-built part of the chain is rotated by a further small randomly chosen angle, in order to 
trigger better sampling. As the last action is nothing else than a pivot move, a full weight 
factor for the backbite move is Wb/w a = (fia sin 0b)/{r sin 9 a ). A single MC cycle contains at 
least N move attempts, randomly selected. A slig htly different version of the backbite move 



has been successfully applied in Refs. |ll|, |l2j, |28| . 



B. Canonical approach: replica exchange 



The parallel tempering (or replica exchange) technique [29|, |30[ is a powerful method 
for sampling in systems with rugged energy landscape. It allows the system to rapidly 
equilibrate and artificially cross energy barriers at low temperatures. Furthermore, the 
method can be easily implemented on a parallel computer. Parallel tempering technique can 
be used with both Monte Carlo and Molecular Dynamics simulations, but in this study, we 
apply this technique to Monte Carlo simulation. The method entails monitoring M canonical 
simulations in parallel at M different temperatures, 7$, i = 1, 2, . . . , M . Each simulation 
corresponds to a replica, or a copy of the system in thermal equilibrium. In individual Monte 
Carlo simulation, new moves are accepted with standard acceptance probabilities given by 



the Metropolis method: 



ri->r. 



min ( 1 , exp 



k B Ti 



(3.2) 



where E^ and E i are the energies of the present and the new conformations, I\ and F i , 
respectively. The replica exchange technique allows the replicas at different temperatures 
to swap with each other without affecting the equilibrium condition at each temperature. 
Specifically, for two replicas, Tj being at Tj and Tj being at Tj, the swap move leads to a 
new state, in which r» is at Tj and Tj is at Tj. The acceptance probability of such a move 
can be derived based on the condition of detailed balance and is given by: 



^(r i ,TO(r J ,T i )^(r J ,T i )cr i ,T i ) = min l,exp 



1 



1 



k B Ti k B Tj 



(Ei - Ej 



(3.3) 



3l| 



The choice of replicas to perform an exchange can be arbitrary, but for a pair of temperatures, 
for which replicas are exchanged, the number of swap move trails must be large enough to 
warrant the statistics. The efficiency of a parallel tempering scheme depends on the number 
of replicas, the set of temperatures to run the simulations, how frequent the swap moves 
are attempted, and is still a matter of debate. It has been suggested that for the best 
performance, the acceptance rate of swap moves must be about 20% 

In our parallel tempering scheme we consider 20 replicas and the temperatures are chosen 
such that they decrease geometrically: T i+1 = aTj, where a = 0.8. The highest reduced 
temperature is k B T\/e = 10, at which the polymer is well poised in the swollen phase. We 
allow replica exchange only between neighboring temperatures and for each replica a swap 
move is attempted every 50 Monte Carlo steps. Standard pivot and crank-shaft move sets 
are used in Monte Carlo simulations. A typical length of the simulations is 10 9 steps per 
replica. 

Results from parallel tempering simulations are equilibrium data and are convenient to 



32] . The latter allows one 



be analyzed using the weighted histogram analysis method 
to estimate the density of states as well as to calculate the thermodynamic averages from 
simulation data at various equilibrium conditions in such a way that minimizes the statistical 
errors where the histograms overlap. 
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IV. RESULTS 



The aim of the calculation is the computation of the density of states (DOS) g(E) as 
remarked. In WL method, g(E) is constructed iteratively, with smaller scale refinements 
made at each level of iteration, controlled by the flatness of energy histogram. We typically 
consider an iteration to reach convergence after 26 — 30 levels of iteration, corresponding 
to a multiplicative factor values of / = 1 — 10~ 8 4- 1 — 10 -9 . This choice is neither unique, 
nor universally accepted, and as this point is crucial for our analysis, it is discussed in some 



details 
Bhatt 



ow. Seaton and coworkers 



16] argued 20 iterations to suffice, while Zhou and 



have additionally shown, that the statistical error of the WL method scales with 
/ as y/1 — f and thus it is of order of 10~ 3 after 20 and 10~ 4 after 26 iteration steps. We 
have explicitly checked this point in our simulations (see Figj2]), where the DOS for a chain 
of N = 16 is reported for both 20 and 30 iterations, with virtually indistinguishable results. 
As a result, the value 26 was used for most of the subsequent simulations in order to keep 
the computational effort low. However, for longest chain length considered in the present 
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FIG. 2. DOS for N = 16, A = 1.5 and E min = -49 after 20 and 30 iterations. The two results 
are practically indistinguishable with a normalized root-mean-square-deviation (NRMSD) between 
the data of 0.002110. 

work (N = 32), we do observe a slight dependence on number of iterations. On the other 
hand, with extending the number of iterations the better quality of data is not guaranteed, 
since the error saturation plays an increasingly important role. This point has been raised 



by several groups 3J, |35| with the rule-of-thumb result that an increase of the number of 
iterations does not necessarily solve the problem since error saturation is an intrinsic feature 



of flatness-controlled WL simulation, 
appears to improve the situation 
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ime-controlled iteration, offered by same authors, 
35|, but it gives disappointing results when applied of 
to a simple hydrophobic-polar (HP) model of protein folding, with some of the low energy 
states resulting unaccessible for a long time, thus providing unsampled regions of DOS. 



Additional recipes were offered by Swetnam and Allen 35|, based on the works of Zhou 
and Su 36j, but we find the original flatness-controlled algorithm to be more reliable, in 



the sense that if the algorithm does not converge within reasonable amount of time, as a 
results of error saturation and/or poor sampling of low-energy states, the correct sampling 
is clearly affected and the run should be discarded. 

An additional crucial step in WL algorithm hinges in the selection of ground state energy. 
As this must be defined at the outset, and is known to drastically affect the low-energy 



behavior of the system 



15 



16], care must be exercised in its selection. At the present time, 



however, there is no universally accepted procedure for off-lattice Wang-Landau simulations, 
and in the present paper we will be following the procedure suggested in Refs. 

HQ 

, that 

has been reported to be reliable in most of the cases. A preliminary run with no low-energy 
cutoff is carried out for a sufficient number of MC steps (10 8 A r in our case). This provides 
an estimate of the minimal ground state energy. In order to avoid poor sampling and large 
computational time, this value is increased of few percents (about 2% in our case), and the 
result is used as the "practical" estimate of the ground state energy (see detailed description 



in Section II of Ref. 



Ill)- 



We have explicitly performed this procedure for chains ranging from N = 4 to N = 128. 
This is depicted in Fig. [31 where the reduced minimum energy per monomer —E min /(Ne) is 
plotted against 1/N. There is a clear trend to saturate toward an estimated —E min / (Ne) ~ 5 
that appears at longer chain lengths. Whether this is related to close packing effects and 
the symmetry of the ground state conformation, is still unclear. Anyway, the lowest energy 
per unit of chain length appears to saturate to a well defined value at longer chain length. 

To illustrate the danger of using an incorrect value to predict the DOS, and quantify 
its dependence with increasing length of the polymer (that is with increasing N), in Fig. H] 
we report the calculation of the heat capacity given by Eq.f l2.9p for two different choices 
of E min /e and different number of monomers N, as a function of the reduced temperature 
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FIG. 3. Plot of the reduced minimum energy per monomer —E m i n /{Ne) as a function of 1/N for 
A = 1.5. The extrapolated value to N — > oo is —E m i n /(Ne) = 4.85. 

ksT/e. Clearly, the low temperatures region is significantly affected by a different choice of 
the ground state, with an error gradually decreasing with N. For instance, when N = 8 cases, 
the difference between E min /e = — 16 and E min /e = — 17 is larger than with E min /e = — 115 
and E min /e = — 117 when N = 32. 
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FIG. 4. Plot of the reduced heat capacity per monomer C{T) / [Nks) as a function of ksT/e for 
three different values (from left to right) of polymer length N = 8, 16 and 32, denoted as a), b) 
and c), correspondingly, and for different ground state energies. In all cases A = 1.5. 



As a prelimi nary step, we have tested our code against exact analytical results valid for 



small N. 



Q,Q, r 



26j, and against previous results using other techniques 



In all cases, we 
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found very good agreement as detailed below. 



Taylor [9 



10], and later Magee et al [26], have computed the DOS analytically for short 



chain lengths N = 4, 5, 6. Using the MATHEMATICA system Ref. 37|, we have reproduced 
the results for tetramers and pentamers at different values of A reported in Tables 1 and 2 
of Ref. 26[, as well as in Figs. 2,3 and 4 of Ref. 10[. With this being done, we have 



then compared results from our simulation code for the same A = 4, 5, always finding an 
excellent agreement. As example of this with N = 5 (pentamer) and A = 1.5, is reported in 
Fig. [5], both for the reduced heat capacity per monomer C{T) / {Nks) and for the internal 
energy per monomer U(T)/(Ne). 
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FIG. 5. Plot of the reduced heat capacity per monomer C{T) / [Nks) and of the internal energy 
per monomer U(T)/N as a function of ksT/e from our WL code and the exact results in the case 
of N = 5. Here A = 1.4, and E m i n /e = —6. 



Further support to the correctness of our WL code stems from a comparison with the 
results by Zhou and Karplus js], who studied the same model using discontinuous molecular 
dynamics (DMD) and canonical MC simulations. Some test runs for small chains N = 4 — 16, 
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are reported in Fig. |6l and can be checked against Fig.3 in Ref. [8[. In all cases, a very good 
agreement is found for the heat capacity per monomer C(T) /(Nks), that is known to be a 
very sensible probe. Note that, in addition to the well-expressed maximum, last three chain 
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FIG. 6. Reduced heat capacity per monomer C(T) / '(Nks) for N = 4,5,6,8, 12, 16 as a function of 
/cgT/e, at A = 1.5. These findings are in very good agreement with Zhou and Karplus (see Fig.3 
in Ref. 



lengths (N = 8,12,16) indicate the presence of two less distinct maxima/plateaux. It is 
worth stressing that in the case N = 16, our calculation is able to probe lower temperature 
regions than in Ref. 8[ , thus enlightening the appearance of the maximum that results blurred 
in Ref. js]. This is because, results from Ref.js] are affected by very large errors at low 
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temperatures (as commonly found in canonical calculations) even for moderate chain lengths. 
For instance, their case N = 16 shows the onset of a large error below ksT/e = 0.3, so lower 
temperatures are in fact cut out (see again their Fig. 3). Conversely, in WL calculation, 
low and high temperatures are equally well sampled and, because of this, we have managed 
to highlight the appearance of peak in the heat capacity at low temperatures, that is an 
indicator of a possible structural transition (see Fig. [6]). WL a ppr oach is therefore extremely 
useful in this respect. As we have mentioned, Zhou and Bhatt 33j estimated the convergence 
error in WL algorithm to scale as a/1 — J, so a conservative estimate of the error is of the 
order 10~ 4 . 

Additional insights within the micro-canonical approach can be obtained by computing 
the first (3(E) and second j(E) derivative of the inverse micro-canonical temperature as 
given in Eqs. (12.31) and (12. 4p . as they are known to be good proxies of structural changes, 
with also the possibility to distinguish between first and second order transitions [27 1. In 
particular, the extrema of the second derivative 'y(E) indicate the corresponding transition 
energies, with negative and positive values associated with second and first order transition, 
respectively. Our results are depicted in FigJTJ where both these quantities are computed 
as a function of the reduced energy per monomer E/(Ne). The clear difference between 
the A = 1.05 and the remaining other values of A is a reflection of the existence of a direct 
coil-crystal transition, without passing through an intermediate globular state, that was 
discussed by Taylor et al 11|, 12[ with much longer chains. 
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FIG. 7. Plot of the reduced first (ej3(E)) and second (e 2 j(E)) derivative of the inverse micro- 
canonical temperature, as a function of the reduced energy per monomer E/(Ne). Considered 
values of A are 1.05, 1.10, 1.20, 1.30. 
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In line with other previous results, we note that the radius of gyration is not a good probe 
for a detailed description of structural changes occurring in the chain at low temperatures, 
since globule - crystal transition does not involve any significant change of average size. This 
can be seen, for instance, in Fig. |H1 where the mean square radius of gyration per monomer 
(Rg) /(Nl 2 ) is plotted as a function of the reduced temperature ksT/e for the same chain 
lengths as above. While it can be clearly seen that each polymer experience a significant 
shrinks in the region of temperatures between 1.0 and 3.0, that is in the same interval where 
a very weak high-temperature peak/plateau is observed on heat capacity curve (see Fig. [6]), 
no noticeable changes occur on further cooling down, as opposed to the heat capacity, that 
shows additional, well pronounced peaks at low temperatures. 

As anticipated, one of our main aims was a careful comparison between micro-canonical 
and canonical approach, in order to assess the strengths and weaknesses of each method in 
the respective domains, and on the reliabilities of results obtained for short chains. Therefore, 
we have cross-checked our results with specialized Monte Carlo simulations in the canonical 
ensemble, using parallel tempering and replica exchange improvements for chain lengths up 
to N = 32. For chain lengths N = 5, 8 and 12 the comparison with exact and WL solutions 
are indistinguishable. The main advantage of the canonical method, as compared to the 
WL counterpart, lays on the fact that we do not have to guess the ground state from the 
outset, as it will naturally emerge as an equilibrium state at sufficiently low temperatures. 
The drawback for longer chains is, of course, that at low temperatures the system becomes 
more and more compact, and a correct sampling becomes increasingly difficult. 

In this parallel calculations, we find a significant difference between the ground state ener- 
gies as computed from the canonical and the Wang-Landau method, at all considered values 
of A. We further note that for temperatures just above the ground state, the curvature of 
temperature dependence of the internal energy in low-temperature region is slightly different 
in the two ensembles, although the absolute values of the energy look quite similar. This 
is reported in Fig. [9] (b) and (c), with representative snapshots of the initial and final state 
stemming from the WL calculations depicted in FigJTUl This difference is magnified in the 
computation of the heat capacity, as shown in FigJTTT where one can clearly see that for the 
cases A = 1.2 (b) and A = 1.3 (c), peak locations and heights differ. We interpret the first 
peak (high temperature) in the heat capacity to correspond to the coil-globule transition, 
and this is supported by the previous results on the mean radius of gyration in Fig. |HJ The 

15 




0.1 1 10 

FIG. 8. Mean square radius of gyration {Rg) /(Nl 2 ) as a function of ksT/e and different length 
chains, and for chain lengths N = 4, 5, 6, 8, 12, 16 at A = 1.5. One can compare these results against 
Fig.l of ij. 

additional peaks appearing at lower temperatures are indications of further globule-globule 
structural changes, not mirrored by the radius of gyration, as remarked. 

It can be easily checked, by matching the energies U(T) corresponding to the temperature 
transitions in Fig. |9l that they nicely match the transition energies reported in FigJT] from 
the micro-canonical approach. 

We note the small discrepancy between the low energy behaviors as obtained from WL and 
from canonical approaches. This is of course always possible in such both have advantages 
and disadvantages that are somewhat complementary one-another. The canonical ensemble 
calculation is more efficient in predicting the correct absolute value of the ground state 
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FIG. 9. Plot of of the internal energy per particle U(T)/(Ne), as a function of ksT/e, for three 
different values of the parameter A: (a) A = 1.05, (b) A = 1.2, (c) A = 1.3. In all cases N = 32. 
Both the WL and the canonical results are reported. 



energy, as the system is naturally driven toward the absolute minimum by the annealing 
process, unlike the WL scheme where this is approximately estimated during the initial 
run. When complemented by replica exchange techniques, allowing a constant swapping of 
conformation between high and low temperatures, the canonical scheme has proven very 
reliable in the correct sampling of configuration space at all temperatures. Of course, the 
sampling becomes increasingly difficult at low temperatures due to the non-swapping moves, 
as remarked. In WL scheme, on the other hand, sampling is equally achieved at all energies 
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E/e=75 E/8=92 E/£=98 

FIG. 10. Representative snapshots of coil (top) and globule (bottom) are shown, for three different 
values of the parameter A: (a) A = 1.05, (b) A = 1.2, (c) A = 1.3. In all cases iV = 32. Energies 
for each of the structures shown on the figure. 



belonging to the chosen interval, including those typically occurring at low temperatures. On 
the other hand, if the lowest energies are not correctly accounted for, one loses an important 
contribution from those states in thermodynamical averages, which can be very well more 
significant at low temperatures. 

Although results obtained with WL and canonical methods slightly differ from one an- 
other at low temperatures, they both provide the same physical picture indicating the ap- 
pearance of two transitions on the phase diagram: coil - molten globule at higher temper- 
atures and molten globule - globule at lower temperatures. In spite of the much shorter 



polymer lengths, our results are in qualitative agreement with those obtained in Refs. |lll.ll2| 
for much longer polymers (N = 128). As in that case, indeed, at values of A close to unit 
unity, there is only one, coil - crystal transition, as indicated by the inset in Fig{TT] (a). 

We have also contrasted results from WL and canonical approach for the average radius 
of gyration per monomer (R 2 ) / '(Nl 2 ) as a function of temperature. This is reported in 
FigH2] for the same parameters as above. Two points are noteworthy. First the critical 
temperature of the coil-globule transition is in agreement with those predicted in Figs j9] and 
ITTI in both approaches. Second, both results indicate that the radius of gyration does not 
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FIG. 11. Plot of the temperature dependence of the reduced heat capacity per particle 
C(T) I '(Nks), for three different values of the parameter A: (a) A = 1.05, (b) A = 1.2, (c) A = 1.3. 
In all cases ./V = 32. Both WL and canonical results are reported. This inset in (a) is a blow-up of 
the most representative part of the figure. 

display any significant change in size beside the first coil-globule transition, as one should 
expect as an indication of globule - crystal transition, again in contrast with the previous 
picture hinging upon the behavior of heat capacity. In this respect, we thus confirm that 
the radius of gyration is not a very good probe of these type of transitions. However, for the 
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FIG. 12. Plot of the temperature dependence of the radius of gyration per particle (R 2 ) /(Nl 2 ), for 
three different values of the parameter A (left-to-right): (a) A = 1.05, (b) A = 1.2, (c) A = 1.3. In 
all cases N = 32. Again, both results from the WL and the canonical approaches are displayed. 



case of A = 1.05 one can notice a sharp change of the radius of gyration near the transition 
temperature T = 0.34, in contrast to the gradual changes of R g in the two other cases. This 
is a manifestation of the first-order like direct transition from coil to compact crystal-like 
phase observed for this value of A. 



V. CONCLUSIONS 

In this paper we have studied the equilibrium statistics of a homopolymer formed by a 
sequence of tangent identical monomers represented by impenetrable hard spheres. Non- 
consecutive spheres, interact via a square-well potential thus driving the collapse of the chain 
at sufficiently low temperatures. Both Wang-Landau micro-canonical and replica-exchange 
canonical calculations were performed for polymers up to N = 32 monomers. We have 
then privileged cross-checking between different approaches over extensive simulations of 
very long chains. In this respect, our approach is complementary to those carried out by 
Taylor et al fill , , where much longer chains (up to N = 256 within a single approach) 
were studied. This comparison allows to uncover the pros and cons of each approach for 
short chains where, presumably, an exhaustive comparison can be carried out. Our results 



are in complete agreement with those from Taylor et al [111 . Il2| . in such we also observe 
evidence of a double transition coil-globule at higher temperatures, and globule - crystal 
transitions at lower temperatures, by working at much shorter polymer lengths, and hence 



with a significant less computational effort involved, that in the work by Taylor et al 111. 1 121] . 
is at the edge of present numerical capability. 
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Several possible perspectives can be envisaged as a continuation of the present work. 



allowing consecutive monomers to inter-penetrate, a local stiffness can be enforced |8j, [26], 
thus allowing for other possible transitions, in addition to those reported above, including 
other morphologies such as helices and tori. A similar effect can be obtained by using a 
tubular chain that breaks the spherical symmetry of the present model 
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. Both these 



models differ from the present one, in such they might allow for a coil-helix transition, 
that cannot be obtained by any spherical symmetric potential. They also require higher 
computational effort and hence must be restricted to small chains only. Work along these 
lines are underway and will be reported in a future publication. 
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